/*
 * Tfunc3D.cpp
 *
 *  Created on: Nov 9, 2011
 *      Author: yi
 */

#include "mbuild/Tfunc3D.h"

Tfunc3D::Tfunc3D() {

}

Tfunc3D::~Tfunc3D() {
}

double Tfunc3D::operator()(double x,double y,double z)
{
	return f(x,y,z);
}

double integrate(Tfunc3D &f,double xl,double xh,double yl,double yh,double zl,double zh,int nx,int ny,int nz)
{
	/* Do integration over the box xl-xh,yl-yh,zl-zh
	 *
	 */
	return integrate_x(f,xl,xh,yl,yh,zl,zh,nx,ny,nz);
}

inline double integrate_z(Tfunc3D &f,double x,double y,double zl,double zh,int nz)
{
	double dz=(zh-zl)/nz;
	double sum=0;
	double d1=0.5-sqrt(3.0)/6.0;
	double d2=0.5+sqrt(3.0)/6.0;
	int k;
	double za,zb,fa,fb;
	for(k=0;k<nz;k++)
	{	za=zl+(k+d1)*dz;
		zb=zl+(k+d2)*dz;
		fa=f(x,y,za);
		fb=f(x,y,zb);
		sum+=fa+fb;
	}
	sum*=dz/2;
	return sum;
}
inline double integrate_y(Tfunc3D &f,double x,double yl,double yh,double zl,double zh,int ny,int nz)
{
	double dy=(yh-yl)/ny;
	double d1=0.5-sqrt(3.0)/6.0;
	double d2=0.5+sqrt(3.0)/6.0;
	int k,i;
	double ya,yb,fa,fb;
	double sum=0;
	for(k=0;k<ny;k++)
	{	ya=yl+(k+d1)*dy;
		yb=yl+(k+d2)*dy;
		fa=integrate_z(f,x,ya,zl,zh,nz);
		fb=integrate_z(f,x,yb,zl,zh,nz);
		sum+=fa+fb;
	}
	sum*=dy/2;
	return sum;
}

inline double integrate_x(Tfunc3D &f,double xl,double xh,double yl,double yh,double zl,double zh,int nx,int ny,int nz)
{
	/* Integration using Quadrature rules
	 *
	 */
	double dx=(xh-xl)/nx;
	double d1=0.5-sqrt(3.0)/6.0;
	double d2=0.5+sqrt(3.0)/6.0;
	int k;
	double xa,xb,fa,fb;
	double sum=0;
	for(k=0;k<nx;k++)
	{	xa=xl+(k+d1)*dx;
		xb=xl+(k+d2)*dx;
		fa=integrate_y(f,xa,yl,yh,zl,zh,ny,nz);
		fb=integrate_y(f,xb,yl,yh,zl,zh,ny,nz);
		sum+=fa+fb;
	}
	sum*=dx/2;
	return sum;
}



